Clean Multiple Regression Code
Hayden Vega and An Bui
2025-05-05
1. Set up
#.........................Load Libraries.........................
library(janitor)
library(tidyverse)
library(ggeffects)
library(sjlabelled)
library(lme4)
library(DHARMa)
library(parameters)
library(effects)
library(MuMIn)
library(waffle)
library(showtext)
library(lmerTest)
library(here)
#.........................Load Fonts.........................
font_add_google(name = "Marcellus", family = "marc")
font_add_google(name = "Josefin Sans", family = "josefin")
font_add_google(name = "Lexend", family = "lex")
#...............Assigning a Standard ggplot Theme............
# ggplot theme
theme_set(theme_bw() +
theme(panel.grid = element_blank(),
plot.background = element_blank()))
# abline color
ref_line_col <- "grey"
# model prediction color
model_col <- "darkblue"
#.......Reading in Data for all Seven Traits and Treatments.........
# Each data frame is read in, then specimen_id column is separated to create new
# columns, then we create an individual ID column from date, site, and individual.
# using read_csv and `here` package to get rid of file paths
height_all <- read_csv(here("data", "preliminary_data", "final_height.csv")) |>
# separate specimen_id column by - to create new date, site, individual, and specimen columns
separate_wider_delim(specimen_id,
delim = "-",
names = c("date", "site", "individual", "specimen"),
cols_remove = FALSE) |>
# create an individual id column
unite("individual_id", date:individual, remove = TRUE, sep = "-")
width_all <- read_csv(here("data", "preliminary_data", "final_width.csv")) |>
separate_wider_delim(specimen_id,
delim = "-", names = c("date", "site", "individual", "specimen"),
cols_remove = FALSE) |>
unite("individual_id", date:individual, remove = TRUE, sep = "-")
perimeter_all <- read_csv(here("data", "preliminary_data", "final_perimeter.csv")) |>
separate_wider_delim(specimen_id,
delim = "-",
names = c("date", "site", "individual", "specimen"),
cols_remove = FALSE) |>
unite("individual_id", date:individual, remove = TRUE, sep = "-")
surf_all <- read_csv(here("data", "preliminary_data", "final_surf.csv")) |>
separate_wider_delim(specimen_id,
delim = "-",
names = c("date", "site", "individual", "specimen"),
cols_remove = FALSE) |>
unite("individual_id", date:individual, remove = TRUE, sep = "-")
tdmc_all <- read_csv(here("data", "preliminary_data", "final_tdmc.csv")) |>
separate_wider_delim(specimen_id,
delim = "-",
names = c("date", "site", "individual", "specimen"),
cols_remove = FALSE) |>
unite("individual_id", date:individual, remove = TRUE, sep = "-")
thickness_all <- read_csv(here("data", "preliminary_data", "final_thickness.csv")) |>
separate_wider_delim(specimen_id,
delim = "-",
names = c("date", "site", "individual", "specimen"),
cols_remove = FALSE) |>
unite("individual_id", date:individual, remove = TRUE, sep = "-")
volume_all <- read_csv(here("data", "preliminary_data", "final_volume.csv")) |>
separate_wider_delim(specimen_id,
delim = "-",
names = c("date", "site", "individual", "specimen"),
cols_remove = FALSE) |>
unite("individual_id", date:individual, remove = TRUE, sep = "-")2. Linear mixed effect models
a. initial values predicted by rehydrate values and/or species
This section includes code to fit linear mixed effect models to
answer the question: how well do rehydrated trait values predict initial
trait values? Each trait is predicted by two models: initial values as a
function of rehydrated values with a species
interaction, or intiial values as a function of rehydrated values
without a species interaction - both models have
individual_id included as a random effect. We view model summaries using
summary().
Height
Width
#..........................Width Models.........................
width_reg_ir <- lmer(width_i ~ width_r * species + (1 | individual_id), data = width_all)
summary(width_reg_ir)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: width_i ~ width_r * species + (1 | individual_id)
## Data: width_all
##
## REML criterion at convergence: 16.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8108 -0.3944 -0.1178 0.3648 5.3792
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0.00000 0.0000
## Residual 0.04315 0.2077
## Number of obs: 162, groups: individual_id, 54
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.00229 0.11672 140.00000 0.020 0.9844
## width_r 1.03003 0.03538 140.00000 29.115 <2e-16 ***
## speciesCC 0.02795 0.15960 140.00000 0.175 0.8613
## speciesCO 0.03873 0.22403 140.00000 0.173 0.8630
## speciesCYOS 0.16350 0.16650 140.00000 0.982 0.3278
## speciesDL 0.17353 0.77701 140.00000 0.223 0.8236
## speciesDP 0.03608 0.13583 140.00000 0.266 0.7909
## speciesEGME 0.47356 0.19897 140.00000 2.380 0.0187 *
## speciesGR -0.12681 0.25299 140.00000 -0.501 0.6170
## speciesNAND 0.07054 0.16759 140.00000 0.421 0.6745
## speciesPOLA 0.12897 0.67780 140.00000 0.190 0.8494
## speciesR 0.04267 0.14411 140.00000 0.296 0.7676
## width_r:speciesCC -0.03103 0.03677 140.00000 -0.844 0.4002
## width_r:speciesCO -0.02300 0.06060 140.00000 -0.379 0.7049
## width_r:speciesCYOS -0.03545 0.04123 140.00000 -0.860 0.3914
## width_r:speciesDL -0.18023 0.65790 140.00000 -0.274 0.7845
## width_r:speciesDP -0.01653 0.04568 140.00000 -0.362 0.7180
## width_r:speciesEGME -0.15968 0.06262 140.00000 -2.550 0.0118 *
## width_r:speciesGR 0.07176 0.10347 140.00000 0.694 0.4891
## width_r:speciesNAND -0.03328 0.04590 140.00000 -0.725 0.4695
## width_r:speciesPOLA -0.06468 0.19565 140.00000 -0.331 0.7415
## width_r:speciesR -0.05309 0.03821 140.00000 -1.389 0.1669
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')width_reg_ir_nospecies <- lmer(width_i ~ width_r + (1 | individual_id), data = width_all)
summary(width_reg_ir_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: width_i ~ width_r + (1 | individual_id)
## Data: width_all
##
## REML criterion at convergence: -38.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0894 -0.3833 -0.1711 0.4370 5.3312
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0.00000 0.0000
## Residual 0.04256 0.2063
## Number of obs: 162, groups: individual_id, 54
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.645e-02 2.502e-02 1.600e+02 3.456 0.000703 ***
## width_r 9.939e-01 4.949e-03 1.600e+02 200.819 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## width_r -0.762
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')### Perimeter
perimeter_reg_ir <- lmer(per_i ~ per_r * species + (1 | individual_id), data = perimeter_all)
summary(perimeter_reg_ir)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: per_i ~ per_r * species + (1 | individual_id)
## Data: perimeter_all
##
## REML criterion at convergence: 1298.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4623 -0.1754 -0.0524 0.1051 5.7167
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 95.52 9.774
## Residual 317.46 17.817
## Number of obs: 153, groups: individual_id, 52
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.94377 9.72785 88.63228 1.228 0.2228
## per_r 0.89901 0.06292 128.62071 14.289 < 2e-16 ***
## speciesCC -12.00315 14.14248 102.00347 -0.849 0.3980
## speciesCO -13.67406 16.07381 104.35426 -0.851 0.3969
## speciesCYOS -6.27899 15.71169 113.52599 -0.400 0.6902
## speciesDP -1.15338 13.35059 94.16302 -0.086 0.9313
## speciesEGME -10.09572 14.24611 108.75008 -0.709 0.4800
## speciesGR -7.51742 16.61084 85.79641 -0.453 0.6520
## speciesNAND -29.82349 34.36413 87.86417 -0.868 0.3878
## speciesPOLA 188.18674 27.83589 131.21107 6.761 4.09e-10 ***
## speciesR -10.42931 13.45056 89.53546 -0.775 0.4402
## per_r:speciesCC 0.12080 0.13574 129.10913 0.890 0.3752
## per_r:speciesCO 0.13346 0.13865 132.86649 0.963 0.3375
## per_r:speciesCYOS 0.10090 0.08751 132.78026 1.153 0.2510
## per_r:speciesDP 0.17766 0.07971 126.37640 2.229 0.0276 *
## per_r:speciesEGME 0.06577 0.20012 107.37549 0.329 0.7431
## per_r:speciesGR 0.09556 0.11385 125.76078 0.839 0.4029
## per_r:speciesNAND 0.35021 0.22788 107.56440 1.537 0.1273
## per_r:speciesPOLA -1.04870 0.13015 110.28003 -8.058 1.01e-12 ***
## per_r:speciesR 0.09930 0.10495 131.89298 0.946 0.3458
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1perimeter_reg_ir_nospecies <- lmer(per_i ~ per_r + (1 | individual_id), data = perimeter_all)
summary(perimeter_reg_ir_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: per_i ~ per_r + (1 | individual_id)
## Data: perimeter_all
##
## REML criterion at convergence: 1434.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3939 -0.2103 -0.1346 0.0058 4.3400
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 72.39 8.508
## Residual 627.45 25.049
## Number of obs: 153, groups: individual_id, 52
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.56937 4.08128 67.48997 1.365 0.177
## per_r 1.00652 0.02865 104.53362 35.138 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## per_r -0.814Surface Area
surf_reg_ir <- lmer(surf_i ~ surf_r * species + (1 | individual_id), data = surf_all)
summary(surf_reg_ir)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: surf_i ~ surf_r * species + (1 | individual_id)
## Data: surf_all
##
## REML criterion at convergence: 1188.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1047 -0.1974 0.0000 0.1264 7.5343
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 81.54 9.03
## Residual 117.26 10.83
## Number of obs: 155, groups: individual_id, 52
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.58126 7.81026 107.08809 -0.202 0.8399
## surf_r 1.12908 0.08684 134.94728 13.002 <2e-16 ***
## speciesCC 15.42369 9.77328 92.22910 1.578 0.1180
## speciesCO 0.89677 10.38244 105.62177 0.086 0.9313
## speciesCYOS 3.29739 11.49479 122.93006 0.287 0.7747
## speciesDP 9.75021 9.20383 103.96275 1.059 0.2919
## speciesEGME 2.18597 9.97535 97.72563 0.219 0.8270
## speciesGR -1.55355 15.32926 82.71308 -0.101 0.9195
## speciesNAND -3.82483 16.69643 128.64243 -0.229 0.8192
## speciesPOLA -4.48724 20.82619 132.22363 -0.215 0.8297
## speciesR 0.15977 10.70600 82.79907 0.015 0.9881
## surf_r:speciesCC -0.16124 0.08978 134.80058 -1.796 0.0747 .
## surf_r:speciesCO 0.05879 0.48435 128.38318 0.121 0.9036
## surf_r:speciesCYOS 0.01105 0.13090 120.77518 0.084 0.9329
## surf_r:speciesDP 0.10850 0.09911 134.94484 1.095 0.2755
## surf_r:speciesEGME -0.10908 0.12497 123.28075 -0.873 0.3844
## surf_r:speciesGR 0.01837 0.17066 111.20076 0.108 0.9145
## surf_r:speciesNAND 0.25954 0.34672 117.31254 0.749 0.4556
## surf_r:speciesPOLA -0.04449 0.17205 109.89074 -0.259 0.7964
## surf_r:speciesR -0.02750 0.15965 94.00692 -0.172 0.8636
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1surf_reg_ir_nospecies <- lmer(surf_i ~ surf_r + (1 | individual_id), data = surf_all)
summary(surf_reg_ir_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: surf_i ~ surf_r + (1 | individual_id)
## Data: surf_all
##
## REML criterion at convergence: 1273
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2860 -0.2405 -0.0934 0.1431 7.0725
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 164.2 12.82
## Residual 128.3 11.33
## Number of obs: 155, groups: individual_id, 52
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.15738 2.38648 80.79972 3.418 0.00099 ***
## surf_r 1.03768 0.01904 145.54351 54.505 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## surf_r -0.534Weight
tdmc_reg_ir <- lmer(weight_i ~ weight_r * species + (1 | individual_id), data = tdmc_all)## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescalingsummary(tdmc_reg_ir)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight_i ~ weight_r * species + (1 | individual_id)
## Data: tdmc_all
##
## REML criterion at convergence: 2556.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2441 -0.2221 -0.0195 0.2035 7.5446
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0 0
## Residual 1382191 1176
## Number of obs: 159, groups: individual_id, 55
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -859.71069 389.57350 139.00000 -2.207 0.02897 *
## weight_r 1.80897 0.03927 139.00000 46.061 < 2e-16 ***
## speciesCC 1948.48497 582.33898 139.00000 3.346 0.00105 **
## speciesCO 796.50963 592.86217 139.00000 1.343 0.18130
## speciesCYOS 1548.97927 816.97129 139.00000 1.896 0.06004 .
## speciesDP 826.85381 515.37765 139.00000 1.604 0.11090
## speciesEGME 983.09662 599.62798 139.00000 1.640 0.10337
## speciesGR -290.19721 917.05077 139.00000 -0.316 0.75214
## speciesNAND 898.08911 1658.39242 139.00000 0.542 0.58900
## speciesPOLA 672.30438 1364.35381 139.00000 0.493 0.62296
## speciesR 885.35159 604.43919 139.00000 1.465 0.14525
## weight_r:speciesCC -0.81629 0.05323 139.00000 -15.335 < 2e-16 ***
## weight_r:speciesCO -0.50096 0.55557 139.00000 -0.902 0.36878
## weight_r:speciesCYOS 0.61326 0.34120 139.00000 1.797 0.07445 .
## weight_r:speciesDP 0.74481 0.17310 139.00000 4.303 3.16e-05 ***
## weight_r:speciesEGME -0.30748 0.58338 139.00000 -0.527 0.59898
## weight_r:speciesGR -0.08245 0.09342 139.00000 -0.883 0.37901
## weight_r:speciesNAND -0.27214 1.17642 139.00000 -0.231 0.81740
## weight_r:speciesPOLA -0.35454 0.62942 139.00000 -0.563 0.57415
## weight_r:speciesR -0.44324 0.23046 139.00000 -1.923 0.05649 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')tdmc_reg_ir_nospecies <- lmer(weight_i ~ weight_r + (1 | individual_id), data = tdmc_all)## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescalingsummary(tdmc_reg_ir_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight_i ~ weight_r + (1 | individual_id)
## Data: tdmc_all
##
## REML criterion at convergence: 2907
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0646 -0.1959 -0.1068 0.1539 6.4238
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 2316648 1522
## Residual 3891255 1973
## Number of obs: 159, groups: individual_id, 55
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 832.53896 295.38612 69.63116 2.818 0.00628 **
## weight_r 1.39800 0.04207 153.61803 33.228 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## weight_r -0.459
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingThickness
thick_reg_ir <- lmer(thick_i ~ thick_r * species + (1 | individual_id), data = thickness_all)
summary(thick_reg_ir)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: thick_i ~ thick_r * species + (1 | individual_id)
## Data: thickness_all
##
## REML criterion at convergence: -96.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2260 -0.1919 -0.0432 0.0915 7.2546
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0.04666 0.2160
## Residual 0.01759 0.1326
## Number of obs: 162, groups: individual_id, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.120542 0.114638 73.651818 1.052 0.2965
## thick_r 0.139866 0.320579 102.193603 0.436 0.6635
## speciesCC 0.010940 0.220029 122.641583 0.050 0.9604
## speciesCO 0.363711 0.331424 132.362361 1.097 0.2745
## speciesCYOS 0.457297 0.237785 136.671480 1.923 0.0565 .
## speciesDL 0.049458 0.459266 139.979514 0.108 0.9144
## speciesDP 0.124849 0.222560 127.567063 0.561 0.5758
## speciesEGME 0.001116 0.220785 123.673763 0.005 0.9960
## speciesGR 0.360171 1.022235 108.311223 0.352 0.7253
## speciesNAND 0.059782 0.373055 139.094404 0.160 0.8729
## speciesPOLA 0.045172 1.568954 99.927148 0.029 0.9771
## speciesR -0.102688 0.194661 68.868970 -0.528 0.5995
## thick_r:speciesCC 0.616505 0.480729 108.357418 1.282 0.2024
## thick_r:speciesCO -0.029348 0.664557 101.790886 -0.044 0.9649
## thick_r:speciesCYOS 0.198255 0.886893 97.868959 0.224 0.8236
## thick_r:speciesDL -0.139866 2.330202 94.656268 -0.060 0.9523
## thick_r:speciesDP 0.495902 1.270180 135.929443 0.390 0.6968
## thick_r:speciesEGME 0.281383 1.275901 110.142813 0.221 0.8259
## thick_r:speciesGR 0.673354 0.797153 126.255219 0.845 0.3999
## thick_r:speciesNAND -0.263202 1.800285 94.824969 -0.146 0.8841
## thick_r:speciesPOLA -0.497009 14.474206 94.516655 -0.034 0.9727
## thick_r:speciesR 0.765868 0.655984 126.848214 1.168 0.2452
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1thick_reg_ir_nospecies <- lmer(thick_i ~ thick_r + (1 | individual_id), data = thickness_all)
summary(thick_reg_ir_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: thick_i ~ thick_r + (1 | individual_id)
## Data: thickness_all
##
## REML criterion at convergence: -63.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2443 -0.2138 -0.0588 0.1191 7.2879
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0.05910 0.2431
## Residual 0.01742 0.1320
## Number of obs: 162, groups: individual_id, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.13422 0.04655 56.48883 2.883 0.00556 **
## thick_r 0.92174 0.08979 82.77441 10.266 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## thick_r -0.656Volume
volume_reg_ir <- lmer(vol_i ~ vol_r * species + (1 | individual_id), data = volume_all)
summary(volume_reg_ir)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: vol_i ~ vol_r * species + (1 | individual_id)
## Data: volume_all
##
## REML criterion at convergence: 680.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1113 -0.2169 -0.0116 0.2568 4.1150
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 1.765 1.329
## Residual 3.864 1.966
## Number of obs: 162, groups: individual_id, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.30390 0.88276 67.59939 1.477 0.144
## vol_r 1.43146 0.06500 137.46633 22.022 < 2e-16 ***
## speciesCC 0.09332 1.30665 68.25035 0.071 0.943
## speciesCO -1.05684 1.30904 58.95997 -0.807 0.423
## speciesCYOS 0.72007 1.50989 99.54528 0.477 0.634
## speciesDL 4.19610 5.94891 121.58941 0.705 0.482
## speciesDP -0.04726 1.13820 72.30213 -0.042 0.967
## speciesEGME -0.97200 1.28596 66.89426 -0.756 0.452
## speciesGR -3.05725 2.23710 70.34650 -1.367 0.176
## speciesNAND -1.15642 2.09422 90.80367 -0.552 0.582
## speciesPOLA -1.11325 2.86861 122.30784 -0.388 0.699
## speciesR -0.78440 1.38258 68.02105 -0.567 0.572
## vol_r:speciesCC -0.52246 0.09625 138.33744 -5.428 2.48e-07 ***
## vol_r:speciesCO -0.69817 1.14178 81.96921 -0.611 0.543
## vol_r:speciesCYOS 0.17432 0.51313 120.21491 0.340 0.735
## vol_r:speciesDL -1.93146 2.40841 100.95555 -0.802 0.424
## vol_r:speciesDP 0.37187 0.28946 139.31565 1.285 0.201
## vol_r:speciesEGME -0.55641 0.79461 123.35933 -0.700 0.485
## vol_r:speciesGR -0.02982 0.20193 105.54459 -0.148 0.883
## vol_r:speciesNAND 0.69625 1.54659 139.62288 0.450 0.653
## vol_r:speciesPOLA -0.02907 1.19204 101.05446 -0.024 0.981
## vol_r:speciesR -0.28390 0.51591 107.55214 -0.550 0.583
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1volume_reg_ir_nospecies <- lmer(vol_i ~ vol_r + (1 | individual_id), data = volume_all)
summary(volume_reg_ir_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: vol_i ~ vol_r + (1 | individual_id)
## Data: volume_all
##
## REML criterion at convergence: 770.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9613 -0.3012 -0.1227 0.2502 5.3131
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 2.812 1.677
## Residual 4.812 2.194
## Number of obs: 162, groups: individual_id, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.15637 0.32749 65.60806 3.531 0.000763 ***
## vol_r 1.21440 0.04462 155.39471 27.216 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## vol_r -0.459b. Initial values predicted by herbarium values
This section includes code to fit linear mixed effect models to
answer the question: how well do herbarium trait values
predict initial trait values? Each trait is predicted by two models:
initial values as a function of herbarium values with a
species interaction and intiial values as a function of herbarium values
without a species interaction - both models have
individual_id included as a random effect. We view model summaries using
summary().
Height
height_reg_ih <- lmer(height_i ~ height_h * species + (1 | individual_id), data = height_all)
summary(height_reg_ih)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: height_i ~ height_h * species + (1 | individual_id)
## Data: height_all
##
## REML criterion at convergence: 752
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3839 -0.3069 -0.0696 0.1275 7.4336
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0.04431 0.2105
## Residual 5.53471 2.3526
## Number of obs: 167, groups: individual_id, 54
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.53223 1.52878 101.18822 0.348 0.7285
## height_h 1.08408 0.10738 117.09633 10.095 <2e-16 ***
## speciesCC -2.62401 2.44318 134.06144 -1.074 0.2847
## speciesCO -0.61692 2.23952 74.89845 -0.275 0.7837
## speciesCYOS 1.36457 3.52820 137.61217 0.387 0.6995
## speciesDL 5.36919 3.03593 144.86101 1.769 0.0791 .
## speciesDP 0.78055 1.92977 96.65881 0.404 0.6868
## speciesEGME -0.31909 1.83980 96.77177 -0.173 0.8627
## speciesGR 1.50188 5.33634 112.08412 0.281 0.7789
## speciesNAND 0.04037 2.76399 141.33926 0.015 0.9884
## speciesPOLA -0.15879 3.19891 143.87363 -0.050 0.9605
## speciesR -0.20268 2.04852 91.93949 -0.099 0.9214
## height_h:speciesCC 0.25621 0.15139 133.99473 1.692 0.0929 .
## height_h:speciesCO -0.01880 0.21540 72.49721 -0.087 0.9307
## height_h:speciesCYOS 0.05317 0.16504 136.64447 0.322 0.7478
## height_h:speciesDL -0.19227 0.15126 144.71505 -1.271 0.2057
## height_h:speciesDP 0.04072 0.13133 110.99089 0.310 0.7571
## height_h:speciesEGME -0.01431 0.12181 126.89901 -0.117 0.9067
## height_h:speciesGR -0.11581 0.25850 104.85633 -0.448 0.6551
## height_h:speciesNAND -0.03487 0.24531 135.14081 -0.142 0.8872
## height_h:speciesPOLA 0.02274 0.25690 120.08051 0.089 0.9296
## height_h:speciesR 0.01836 0.16646 93.30354 0.110 0.9124
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#This model describes height_i as a function of height_h, with a species interaction term and random effect of individual_id
height_reg_ih_nospecies <- lmer(height_i ~ height_h + (1 | individual_id), data = height_all)
summary(height_reg_ih_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: height_i ~ height_h + (1 | individual_id)
## Data: height_all
##
## REML criterion at convergence: 783.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9741 -0.4023 -0.1279 0.1446 7.3759
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0.8207 0.9059
## Residual 5.5105 2.3474
## Number of obs: 167, groups: individual_id, 54
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.25808 0.46707 115.71914 0.553 0.582
## height_h 1.13738 0.02731 139.85308 41.641 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## height_h -0.880#This model describes height_i as a function of height_h with a random effect of individual_id and without a species interaction termWidth
width_reg_ih <- lmer(width_i ~ width_h * species + (1 | individual_id), data = width_all)
summary(width_reg_ih)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: width_i ~ width_h * species + (1 | individual_id)
## Data: width_all
##
## REML criterion at convergence: 220.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4570 -0.3541 -0.0043 0.4800 3.7396
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0.02578 0.1606
## Residual 0.17125 0.4138
## Number of obs: 162, groups: individual_id, 54
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.272115 0.236478 126.740300 1.151 0.252021
## width_h 1.151373 0.083930 138.419124 13.718 < 2e-16 ***
## speciesCC -0.620575 0.342832 123.037547 -1.810 0.072714 .
## speciesCO 0.446118 0.415375 139.326962 1.074 0.284674
## speciesCYOS 0.471786 0.340900 109.731366 1.384 0.169185
## speciesDL 0.246742 1.344039 112.207950 0.184 0.854673
## speciesDP -0.523920 0.296015 118.418792 -1.770 0.079316 .
## speciesEGME 0.284999 0.408095 132.278112 0.698 0.486176
## speciesGR -0.122673 0.502592 136.144809 -0.244 0.807536
## speciesNAND -0.126218 0.356086 102.952743 -0.354 0.723719
## speciesPOLA -0.433035 1.673731 107.929898 -0.259 0.796341
## speciesR -0.123335 0.301637 110.557970 -0.409 0.683415
## width_h:speciesCC 0.182215 0.088658 138.359422 2.055 0.041735 *
## width_h:speciesCO -0.288921 0.122342 134.718497 -2.362 0.019630 *
## width_h:speciesCYOS -0.090180 0.097947 139.835801 -0.921 0.358791
## width_h:speciesDL -0.373692 1.575169 100.246932 -0.237 0.812956
## width_h:speciesDP 0.467741 0.135512 135.940672 3.452 0.000743 ***
## width_h:speciesEGME -0.196477 0.145111 139.328362 -1.354 0.177934
## width_h:speciesGR 0.008776 0.243724 134.944567 0.036 0.971328
## width_h:speciesNAND -0.053593 0.107060 134.929958 -0.501 0.617479
## width_h:speciesPOLA 0.158052 0.600841 101.038027 0.263 0.793046
## width_h:speciesR -0.028081 0.091096 139.271757 -0.308 0.758347
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1width_reg_ih_nospecies <- lmer(width_i ~ width_h + (1 | individual_id), data = width_all)
summary(width_reg_ih_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: width_i ~ width_h + (1 | individual_id)
## Data: width_all
##
## REML criterion at convergence: 281.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7952 -0.4539 0.0449 0.3985 4.8626
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0.1262 0.3552
## Residual 0.2309 0.4805
## Number of obs: 162, groups: individual_id, 54
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.08737 0.08547 98.88700 1.022 0.309
## width_h 1.20652 0.01892 154.81287 63.760 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## width_h -0.686Perimeter
perimeter_reg_ih <- lmer(per_i ~ per_h * species + (1 | individual_id), data = perimeter_all)
summary(perimeter_reg_ih)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: per_i ~ per_h * species + (1 | individual_id)
## Data: perimeter_all
##
## REML criterion at convergence: 1361.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0987 -0.2817 -0.0422 0.1080 4.7679
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 225.6 15.02
## Residual 501.2 22.39
## Number of obs: 153, groups: individual_id, 52
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.59684 13.77863 88.26129 0.188 0.850942
## per_h 1.19298 0.10995 125.45680 10.850 < 2e-16 ***
## speciesCC -6.82156 19.70908 103.74382 -0.346 0.729959
## speciesCO -19.35573 23.91988 114.13986 -0.809 0.420089
## speciesCYOS 6.64731 21.07977 109.76272 0.315 0.753102
## speciesDP 34.33115 17.91760 88.87197 1.916 0.058572 .
## speciesEGME -2.02207 19.61561 105.50238 -0.103 0.918091
## speciesGR 1.05566 22.86537 80.61552 0.046 0.963290
## speciesNAND -19.63704 47.84044 78.89181 -0.410 0.682574
## speciesPOLA 147.31861 41.95239 132.97267 3.512 0.000609 ***
## speciesR -1.44377 18.60890 89.93418 -0.078 0.938330
## per_h:speciesCC 0.11400 0.22636 126.28471 0.504 0.615398
## per_h:speciesCO 0.18647 0.24758 132.38146 0.753 0.452687
## per_h:speciesCYOS -0.11782 0.13820 132.81018 -0.853 0.395437
## per_h:speciesDP 0.10835 0.13834 122.61000 0.783 0.434995
## per_h:speciesEGME -0.09324 0.29545 110.83387 -0.316 0.752910
## per_h:speciesGR 0.03776 0.19031 116.97364 0.198 0.843049
## per_h:speciesNAND 0.20992 0.35881 95.46443 0.585 0.559902
## per_h:speciesPOLA -1.03952 0.28146 106.94185 -3.693 0.000351 ***
## per_h:speciesR 0.09234 0.17917 132.89084 0.515 0.607163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1perimeter_reg_ih_nospecies <- lmer(per_i ~ per_h + (1 | individual_id), data = perimeter_all)
summary(perimeter_reg_ih_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: per_i ~ per_h + (1 | individual_id)
## Data: perimeter_all
##
## REML criterion at convergence: 1472.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8908 -0.3120 -0.0885 0.1870 4.6894
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 563.2 23.73
## Residual 573.3 23.94
## Number of obs: 153, groups: individual_id, 52
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 10.01440 5.75140 103.30764 1.741 0.0846 .
## per_h 1.21717 0.04487 150.79674 27.126 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## per_h -0.740Surface Area
surf_reg_ih <- lmer(surf_i ~ surf_h * species + (1 | individual_id), data = surf_all)
summary(surf_reg_ih)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: surf_i ~ surf_h * species + (1 | individual_id)
## Data: surf_all
##
## REML criterion at convergence: 1146.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0489 -0.2322 0.0118 0.2319 3.3970
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 99.13 9.956
## Residual 80.88 8.993
## Number of obs: 155, groups: individual_id, 52
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.463120 7.163540 101.638932 -0.065 0.94858
## surf_h 1.652726 0.110907 133.491909 14.902 < 2e-16 ***
## speciesCC 1.695839 9.223655 86.168702 0.184 0.85456
## speciesCO 2.335557 9.411117 97.867238 0.248 0.80452
## speciesCYOS 2.814195 10.432608 112.263576 0.270 0.78785
## speciesDP 15.162740 8.453638 97.601006 1.794 0.07597 .
## speciesEGME 1.351731 9.298801 88.467416 0.145 0.88475
## speciesGR -1.999914 14.572854 81.746029 -0.137 0.89118
## speciesNAND 1.678109 14.062657 115.197954 0.119 0.90522
## speciesPOLA -87.262355 25.937553 132.709192 -3.364 0.00100 **
## speciesR 0.702535 10.027528 79.063684 0.070 0.94432
## surf_h:speciesCC -0.042041 0.115345 132.677458 -0.364 0.71608
## surf_h:speciesCO -0.390872 0.504958 112.882983 -0.774 0.44051
## surf_h:speciesCYOS -0.042381 0.160204 119.030005 -0.265 0.79182
## surf_h:speciesDP 0.007247 0.125121 133.330151 0.058 0.95390
## surf_h:speciesEGME -0.435098 0.142602 123.731498 -3.051 0.00279 **
## surf_h:speciesGR -0.212427 0.206861 113.434524 -1.027 0.30665
## surf_h:speciesNAND -0.064651 0.338893 108.339360 -0.191 0.84906
## surf_h:speciesPOLA 2.188812 0.450745 99.174844 4.856 4.47e-06 ***
## surf_h:speciesR -0.245841 0.197231 105.183512 -1.246 0.21536
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1surf_reg_ih_nospecies <- lmer(surf_i ~ surf_h + (1 | individual_id), data = surf_all)
summary(surf_reg_ih_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: surf_i ~ surf_h + (1 | individual_id)
## Data: surf_all
##
## REML criterion at convergence: 1261.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8828 -0.3740 -0.0307 0.3144 3.5874
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 187.6 13.70
## Residual 111.4 10.56
## Number of obs: 155, groups: individual_id, 52
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.33827 2.47842 81.86394 0.943 0.348
## surf_h 1.59703 0.02794 140.02936 57.160 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## surf_h -0.532Weight
tdmc_reg_ih <- lmer(weight_i ~ weight_h * species + (1 | individual_id), data = tdmc_all)## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescalingsummary(tdmc_reg_ih)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight_i ~ weight_h * species + (1 | individual_id)
## Data: tdmc_all
##
## REML criterion at convergence: 2630.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9009 -0.1250 -0.0089 0.1327 7.8408
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 466323 682.9
## Residual 2460787 1568.7
## Number of obs: 159, groups: individual_id, 55
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -112.6370 598.8637 71.1358 -0.188 0.851346
## weight_h 7.8059 0.2498 138.8975 31.246 < 2e-16 ***
## speciesCC 205.9332 927.2423 77.4660 0.222 0.824827
## speciesCO -23.2533 954.1125 74.7643 -0.024 0.980621
## speciesCYOS 227.2907 1276.3765 130.1371 0.178 0.858941
## speciesDP 155.9781 798.6326 70.8440 0.195 0.845712
## speciesEGME 168.0942 940.4295 83.5829 0.179 0.858573
## speciesGR -1545.9601 1472.0010 68.2519 -1.050 0.297311
## speciesNAND 2311.3346 1718.3377 134.9808 1.345 0.180849
## speciesPOLA 129.8183 1807.9914 117.0949 0.072 0.942882
## speciesR 223.2403 914.6747 80.4121 0.244 0.807802
## weight_h:speciesCC -3.1801 0.3491 136.7897 -9.110 9.01e-16 ***
## weight_h:speciesCO -5.3849 1.5817 117.6456 -3.405 0.000907 ***
## weight_h:speciesCYOS -1.3612 1.2933 138.7170 -1.053 0.294373
## weight_h:speciesDP -1.0085 0.7308 112.6685 -1.380 0.170294
## weight_h:speciesEGME -0.6318 3.8956 133.7570 -0.162 0.871404
## weight_h:speciesGR -3.6201 0.3985 121.2700 -9.085 2.39e-15 ***
## weight_h:speciesNAND -8.2313 3.9984 115.6620 -2.059 0.041773 *
## weight_h:speciesPOLA -0.1031 4.1748 100.1666 -0.025 0.980344
## weight_h:speciesR -3.4064 1.1599 96.2291 -2.937 0.004151 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## fit warnings:
## Some predictor variables are on very different scales: consider rescalingtdmc_reg_ih_nospecies <- lmer(weight_i ~ weight_h + (1 | individual_id), data = tdmc_all)## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescalingsummary(tdmc_reg_ih_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight_i ~ weight_h + (1 | individual_id)
## Data: tdmc_all
##
## REML criterion at convergence: 2959.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8499 -0.1521 -0.0275 0.1360 7.5307
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 4023996 2006
## Residual 5213112 2283
## Number of obs: 159, groups: individual_id, 55
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 157.2287 382.9288 47.7765 0.411 0.683
## weight_h 5.2341 0.1939 139.7789 26.997 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## weight_h -0.504
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingThickness
thick_reg_ih <- lmer(thick_i ~ thick_h * species + (1 | individual_id), data = thickness_all)
summary(thick_reg_ih)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: thick_i ~ thick_h * species + (1 | individual_id)
## Data: thickness_all
##
## REML criterion at convergence: -112.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5330 -0.2257 -0.0360 0.1393 7.2536
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0.04748 0.2179
## Residual 0.01649 0.1284
## Number of obs: 162, groups: individual_id, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.11225 0.11567 76.23562 0.970 0.3349
## thick_h 0.47944 0.86955 107.00347 0.551 0.5825
## speciesCC 0.31780 0.16174 63.87740 1.965 0.0538 .
## speciesCO 0.31043 0.35703 110.60168 0.869 0.3865
## speciesCYOS 0.16492 0.30022 139.77289 0.549 0.5837
## speciesDL 0.05775 0.73919 121.33255 0.078 0.9379
## speciesDP 0.50291 0.22373 130.49152 2.248 0.0263 *
## speciesEGME 0.03790 0.18241 90.42927 0.208 0.8359
## speciesGR 3.77862 0.82284 139.99471 4.592 9.68e-06 ***
## speciesNAND -0.01503 0.92112 134.73544 -0.016 0.9870
## speciesPOLA -0.06745 0.38776 128.17496 -0.174 0.8622
## speciesR -0.10166 0.20280 60.71286 -0.501 0.6180
## thick_h:speciesCC -0.31147 0.90024 106.65416 -0.346 0.7300
## thick_h:speciesCO -0.24292 1.08623 133.32579 -0.224 0.8234
## thick_h:speciesCYOS 1.39921 1.55498 109.97012 0.900 0.3702
## thick_h:speciesDL -0.47944 5.44469 98.06901 -0.088 0.9300
## thick_h:speciesDP -2.79756 1.73680 129.77461 -1.611 0.1097
## thick_h:speciesEGME 0.13574 2.48680 105.53951 0.055 0.9566
## thick_h:speciesGR -2.41704 1.10667 123.25325 -2.184 0.0309 *
## thick_h:speciesNAND 0.28492 11.28077 137.55240 0.025 0.9799
## thick_h:speciesPOLA 1.40537 6.71613 97.98956 0.209 0.8347
## thick_h:speciesR 1.21838 1.43063 115.64928 0.852 0.3962
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1thick_reg_ih_nospecies <- lmer(thick_i ~ thick_h + (1 | individual_id), data = thickness_all)
summary(thick_reg_ih_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: thick_i ~ thick_h + (1 | individual_id)
## Data: thickness_all
##
## REML criterion at convergence: -47.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1179 -0.2023 -0.0820 0.1024 6.7408
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0.05453 0.2335
## Residual 0.02079 0.1442
## Number of obs: 162, groups: individual_id, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.20382 0.04228 47.93657 4.821 1.48e-05 ***
## thick_h 0.95958 0.09883 75.47472 9.710 6.32e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## thick_h -0.588Volume
volume_reg_ih <- lmer(vol_i ~ vol_h * species + (1 | individual_id), data = volume_all)
summary(volume_reg_ih)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: vol_i ~ vol_h * species + (1 | individual_id)
## Data: volume_all
##
## REML criterion at convergence: 722.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8208 -0.2484 -0.0213 0.2188 4.4006
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 0.3959 0.6292
## Residual 7.2424 2.6912
## Number of obs: 162, groups: individual_id, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.2386 1.0413 97.6340 -0.229 0.819210
## vol_h 6.1714 0.3727 139.9985 16.557 < 2e-16 ***
## speciesCC 0.1428 1.5687 96.6246 0.091 0.927664
## speciesCO 0.4539 1.4577 87.4885 0.311 0.756233
## speciesCYOS 0.8583 2.2830 139.9888 0.376 0.707525
## speciesDL 5.2386 4.8174 131.1702 1.087 0.278835
## speciesDP 0.9271 1.3266 87.0716 0.699 0.486484
## speciesEGME 0.7993 1.4005 79.3197 0.571 0.569798
## speciesGR -1.5082 2.7201 107.8951 -0.554 0.580414
## speciesNAND 1.2685 2.5855 133.2270 0.491 0.624496
## speciesPOLA 3.7132 3.7907 139.3476 0.980 0.328996
## speciesR 0.3816 1.7283 81.9842 0.221 0.825814
## vol_h:speciesCC -1.5920 0.5894 128.1774 -2.701 0.007845 **
## vol_h:speciesCO -5.3670 1.3987 124.1825 -3.837 0.000197 ***
## vol_h:speciesCYOS -2.1893 1.7567 117.9363 -1.246 0.215151
## vol_h:speciesDL -6.6714 3.3170 106.3239 -2.011 0.046827 *
## vol_h:speciesDP -2.0224 0.8121 108.2204 -2.490 0.014280 *
## vol_h:speciesEGME -4.9200 2.0074 136.3490 -2.451 0.015517 *
## vol_h:speciesGR -3.6149 0.5987 137.4855 -6.038 1.38e-08 ***
## vol_h:speciesNAND -3.9645 4.1135 127.1705 -0.964 0.336993
## vol_h:speciesPOLA -7.6968 7.0171 105.8129 -1.097 0.275191
## vol_h:speciesR -2.9773 1.8410 68.9227 -1.617 0.110406
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1volume_reg_ih_nospecies <- lmer(vol_i ~ vol_h + (1 | individual_id), data = volume_all)
summary(volume_reg_ih_nospecies)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: vol_i ~ vol_h + (1 | individual_id)
## Data: volume_all
##
## REML criterion at convergence: 888
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0468 -0.2718 -0.0058 0.2196 5.9598
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_id (Intercept) 6.823 2.612
## Residual 9.825 3.135
## Number of obs: 162, groups: individual_id, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.2726 0.5367 57.1506 0.508 0.613
## vol_h 3.7867 0.2330 142.9514 16.252 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## vol_h -0.573#3. Comparing models using AIC
This section includes code comparing the AIC values for models with a species interaction term and without a species interaction term to answer the question: Do models with a species interaction improve the model’s fit? All trait models for rehydrated and herbarium traits are tested, and the results are commented below each model.
a. Comparisons for models with rehydrated traits as predictors
Height
#..........................Comparisons for Height Models........................
model.sel(height_reg_ir, height_reg_ir_nospecies)## Model selection table
## (Int) hgh_r spc hgh_r:spc df logLik AICc delta
## height_reg_ir_nospecies 0.04065 1.024 4 -169.859 348.0 0.00
## height_reg_ir -0.27440 1.029 + + 24 -169.996 396.6 48.66
## weight
## height_reg_ir_nospecies 1
## height_reg_ir 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#nospecies is betterWidth
#..........................Comparisons for Width Models.........................
model.sel(width_reg_ir, width_reg_ir_nospecies)## Model selection table
## (Int) spc wdt_r spc:wdt_r df logLik AICc delta
## width_reg_ir_nospecies 0.08645 0.9939 4 19.245 -30.2 0.00
## width_reg_ir 0.00229 + 1.0300 + 24 -8.349 73.5 103.69
## weight
## width_reg_ir_nospecies 1
## width_reg_ir 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#nospecies is betterPerimeter
#..........................Comparisons for Perimeter Models......................
model.sel(perimeter_reg_ir, perimeter_reg_ir_nospecies)## Model selection table
## (Int) per_r spc per_r:spc df logLik AICc delta
## perimeter_reg_ir 11.940 0.899 + + 22 -649.166 1350.1 0.00
## perimeter_reg_ir_nospecies 5.569 1.007 4 -717.235 1442.7 92.62
## weight
## perimeter_reg_ir 1
## perimeter_reg_ir_nospecies 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#with species is betterSurface Area
#..........................Comparisons for Surface Area Models..................
model.sel(surf_reg_ir, surf_reg_ir_nospecies)## Model selection table
## (Int) spc srf_r spc:srf_r df logLik AICc delta
## surf_reg_ir -1.581 + 1.129 + 22 -594.245 1240.2 0.00
## surf_reg_ir_nospecies 8.157 1.038 4 -636.504 1281.3 41.12
## weight
## surf_reg_ir 1
## surf_reg_ir_nospecies 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#species is betterWeight
#..........................Comparisons for Weight Models........................
model.sel(tdmc_reg_ir, tdmc_reg_ir_nospecies)## Model selection table
## (Int) spc wgh_r spc:wgh_r df logLik AICc delta
## tdmc_reg_ir -859.7 + 1.809 + 22 -1278.379 2608.2 0.00
## tdmc_reg_ir_nospecies 832.5 1.398 4 -1453.517 2915.3 307.09
## weight
## tdmc_reg_ir 1
## tdmc_reg_ir_nospecies 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#species is betterThickness
#..........................Comparisons for Thickness Models.....................
model.sel(thick_reg_ir, thick_reg_ir_nospecies)## Model selection table
## (Int) spc thc_r spc:thc_r df logLik AICc delta weight
## thick_reg_ir_nospecies 0.1342 0.9217 4 31.946 -55.6 0.00 1
## thick_reg_ir 0.1205 + 0.1399 + 24 48.267 -39.8 15.86 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#nospecies is better, but rehydrated thickness it's not a good predictor overallVolume
#..........................Comparisons for Volume Models........................
model.sel(volume_reg_ir, volume_reg_ir_nospecies)## Model selection table
## (Int) spc vol_r spc:vol_r df logLik AICc delta
## volume_reg_ir 1.304 + 1.431 + 24 -340.348 737.5 0.00
## volume_reg_ir_nospecies 1.156 1.214 4 -385.284 778.8 41.37
## weight
## volume_reg_ir 1
## volume_reg_ir_nospecies 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#with species is better##b. Comparisons for models with Herbarium traits as a predictor variable
Height
#..........................Comparisons for Height Models........................
model.sel(height_reg_ih, height_reg_ih_nospecies)## Model selection table
## (Int) hgh_h spc hgh_h:spc df logLik AICc delta
## height_reg_ih_nospecies 0.2581 1.137 4 -391.826 791.9 0.00
## height_reg_ih 0.5322 1.084 + + 24 -376.005 808.5 16.56
## weight
## height_reg_ih_nospecies 1
## height_reg_ih 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#nospecies is betterWidth
#..........................Comparisons for Width Models........................
model.sel(width_reg_ih, width_reg_ih_nospecies)## Model selection table
## (Int) spc wdt_h spc:wdt_h df logLik AICc delta
## width_reg_ih 0.27210 + 1.151 + 24 -110.139 277.0 0.0
## width_reg_ih_nospecies 0.08737 1.207 4 -140.689 289.6 12.6
## weight
## width_reg_ih 0.998
## width_reg_ih_nospecies 0.002
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#species is betterPerimeter
#.......................Comparisons for Perimeter Models........................
model.sel(perimeter_reg_ih, perimeter_reg_ih_nospecies)## Model selection table
## (Int) per_h spc per_h:spc df logLik AICc delta
## perimeter_reg_ih 2.597 1.193 + + 22 -680.895 1413.6 0.00
## perimeter_reg_ih_nospecies 10.010 1.217 4 -736.285 1480.8 67.27
## weight
## perimeter_reg_ih 1
## perimeter_reg_ih_nospecies 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#with species is betterSurface Area
#....................Comparisons for Surface Area Models........................
model.sel(surf_reg_ih, surf_reg_ih_nospecies)## Model selection table
## (Int) spc srf_h spc:srf_h df logLik AICc delta
## surf_reg_ih -0.4631 + 1.653 + 22 -573.459 1198.6 0.00
## surf_reg_ih_nospecies 2.3380 1.597 4 -630.864 1270.0 71.41
## weight
## surf_reg_ih 1
## surf_reg_ih_nospecies 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#species is betterWeight
#..........................Comparisons for Weight Models........................
model.sel(tdmc_reg_ih, tdmc_reg_ih_nospecies)## Model selection table
## (Int) spc wgh_h spc:wgh_h df logLik AICc delta
## tdmc_reg_ih -112.6 + 7.806 + 22 -1315.397 2682.2 0.00
## tdmc_reg_ih_nospecies 157.2 5.234 4 -1479.581 2967.4 285.19
## weight
## tdmc_reg_ih 1
## tdmc_reg_ih_nospecies 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#species is betterThickness
#.......................Comparisons for Thickness Models........................
model.sel(thick_reg_ih, thick_reg_ih_nospecies)## Model selection table
## (Int) spc thc_h spc:thc_h df logLik AICc delta weight
## thick_reg_ih 0.1123 + 0.4794 + 24 56.121 -55.5 0.00 1
## thick_reg_ih_nospecies 0.2038 0.9596 4 23.776 -39.3 16.19 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#species is better, but doesn't really matter cause it's not a good predictorVolume
#..........................Comparisons for Volume Models........................
model.sel(volume_reg_ih, volume_reg_ih_nospecies)## Model selection table
## (Int) spc vol_h spc:vol_h df logLik AICc delta
## volume_reg_ih -0.2386 + 6.171 + 24 -361.033 778.8 0.00
## volume_reg_ih_nospecies 0.2726 3.787 4 -443.999 896.3 117.43
## weight
## volume_reg_ih 1
## volume_reg_ih_nospecies 0
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | individual_id#with species is better4. Calculating R squared values for the best fitting models
This section includes code that calculates the marginal and conditional R squared values of the best fitting models from the section above. We calculate the R squared value using MuMIn::r.squaredGLMM(model_name) from the MuMIn package.
##a. Calculating R squared values for models with rehydrated traits as a predictor
Height
#..................R^2 for Height Model without species interaction term........
MuMIn::r.squaredGLMM(height_reg_ir_nospecies)## R2m R2c
## [1,] 0.9933504 0.9959516Width
#..................R^2 for Width Model without species interaction term........
MuMIn::r.squaredGLMM(width_reg_ir_nospecies)## R2m R2c
## [1,] 0.9960237 0.9960237Perimeter
#...............R^2 for Perimeter Model with species interaction term........
MuMIn::r.squaredGLMM(perimeter_reg_ir)## R2m R2c
## [1,] 0.9417458 0.9552199Surface Area
#............R^2 for Surface Area Model with species interaction term........
MuMIn::r.squaredGLMM(surf_reg_ir)## R2m R2c
## [1,] 0.9630645 0.9782143Weight
#..................R^2 for Weight Model with species interaction term........
MuMIn::r.squaredGLMM(tdmc_reg_ir)## R2m R2c
## [1,] 0.9728602 0.9728602Thickness
#...............R^2 for Thickness Model without species interaction term........
MuMIn::r.squaredGLMM(thick_reg_ir_nospecies)## R2m R2c
## [1,] 0.4802988 0.8816962Volume
#..................R^2 for Volume Model with species interaction term........
MuMIn::r.squaredGLMM(volume_reg_ir)## R2m R2c
## [1,] 0.8780037 0.9162553##b. Calculating R squared values for models with Herbarium traits as predictor variables
Height
#..................R^2 for Height Model without species interaction term........
MuMIn::r.squaredGLMM(height_reg_ih_nospecies)## R2m R2c
## [1,] 0.9182208 0.9288213Width
#..................R^2 for Width Model with species interaction term............
MuMIn::r.squaredGLMM(width_reg_ih)## R2m R2c
## [1,] 0.9814827 0.9839057Perimeter
#...............R^2 for Perimeter Model with species interaction term........
MuMIn::r.squaredGLMM(perimeter_reg_ih)## R2m R2c
## [1,] 0.8997358 0.9308613Surface Area
#............R^2 for Surface Area Model with species interaction term........
MuMIn::r.squaredGLMM(surf_reg_ih)## R2m R2c
## [1,] 0.9662619 0.9848414Weight
#..................R^2 for Weight Model with species interaction term........
MuMIn::r.squaredGLMM(tdmc_reg_ih)## R2m R2c
## [1,] 0.9435099 0.9525094Thickness
#..................R^2 for Thickness Model with species interaction term........
MuMIn::r.squaredGLMM(thick_reg_ih)## R2m R2c
## [1,] 0.6413815 0.9075627Volume
#..................R^2 for Volume Model with species interaction term...........
MuMIn::r.squaredGLMM(volume_reg_ih)## R2m R2c
## [1,] 0.8312203 0.8399674#5. Saving models as predictions for later plotting
This section includes code to save the model predictions of the best fitting models as data frames using ggpredict(). This includes predictions for all traits and both rehydrate and herbarium treatments.
a. Saving predictions for models with rehydrated traits as predictor variables
Height
#........Predictions for the Height Model Without Species........
height_ir_predictions <- ggpredict(height_reg_ir_nospecies, terms = "height_r[0:45]") |>
#saving predictions for height_r values from 0-45
mutate(trait = "height") #adding a column to label the type of trait values
height_ir_predictions## # Predicted values of height_i
##
## height_r | Predicted | 95% CI | trait
## -------------------------------------------
## 0 | 0.04 | -0.23, 0.31 | height
## 1 | 1.06 | 0.81, 1.32 | height
## 2 | 2.09 | 1.84, 2.34 | height
## 3 | 3.11 | 2.87, 3.35 | height
## 4 | 4.14 | 3.91, 4.36 | height
## 5 | 5.16 | 4.94, 5.38 | height
## 6 | 6.18 | 5.98, 6.39 | height
##
## Adjusted for:
## * individual_id = 0 (population-level)Width
#........Predictions for the Width Model Without Species........
width_ir_predictions <- ggpredict(width_reg_ir_nospecies, terms = "width_r[0:30]") |> #max width = 24
mutate(trait = "width")
# look at the data frame
width_ir_predictions## # Predicted values of width_i
##
## width_r | Predicted | 95% CI | trait
## ----------------------------------------
## 0 | 0.09 | 0.04, 0.14 | width
## 1 | 1.08 | 1.04, 1.12 | width
## 2 | 2.07 | 2.04, 2.11 | width
## 3 | 3.07 | 3.04, 3.10 | width
## 4 | 4.06 | 4.03, 4.09 | width
## 5 | 5.06 | 5.02, 5.09 | width
## 6 | 6.05 | 6.01, 6.09 | width
##
## Adjusted for:
## * individual_id = 0 (population-level)Perimeter
#........Predictions for the Height Model With Species........
perimeter_ir_predictions <- ggpredict(perimeter_reg_ir, terms = "per_r[0:375]") |>
mutate(trait = "perimeter")
#...........Incorporating Species as a Prediction Term...........
perimeter_ir_predictions_species <- ggpredict(perimeter_reg_ir, terms = c("per_r", "species")) |>
as_tibble() |> #converting to a tibble
rename(species = group)
# look at the data frame
perimeter_ir_predictions## # Predicted values of per_i
##
## per_r | Predicted | 95% CI | trait
## --------------------------------------------
## 0 | 11.94 | -7.30, 31.19 | perimeter
## 1 | 12.84 | -6.30, 31.99 | perimeter
## 2 | 13.74 | -5.30, 32.79 | perimeter
## 3 | 14.64 | -4.30, 33.59 | perimeter
## 4 | 15.54 | -3.31, 34.39 | perimeter
## 5 | 16.44 | -2.31, 35.19 | perimeter
## 6 | 17.34 | -1.31, 35.99 | perimeter
##
## Adjusted for:
## * species = BF
## * individual_id = 0 (population-level)Surface Area
#........Predictions for the Height Model With Species........
surf_ir_predictions <- ggpredict(surf_reg_ir, terms = "surf_r[0:430]") |>
mutate(trait = "surface area")
#...........Incorporating Species as a Prediction Term...........
surf_ir_predictions_species <- ggpredict(surf_reg_ir, terms = c("surf_r", "species")) |>
as_tibble() |>
rename(species = group)
# look at the data frame
surf_ir_predictions## # Predicted values of surf_i
##
## surf_r | Predicted | 95% CI | trait
## -------------------------------------------------
## 0 | -1.58 | -17.03, 13.87 | surface area
## 1 | -0.45 | -15.76, 14.86 | surface area
## 2 | 0.68 | -14.49, 15.85 | surface area
## 3 | 1.81 | -13.22, 16.84 | surface area
## 4 | 2.94 | -11.96, 17.83 | surface area
## 5 | 4.06 | -10.69, 18.82 | surface area
## 6 | 5.19 | -9.42, 19.81 | surface area
##
## Adjusted for:
## * species = BF
## * individual_id = 0 (population-level)Weight
#........Predictions for the Weight Model With Species........
tdmc_ir_predictions <- ggpredict(tdmc_reg_ir, terms = "weight_r[0:29050]") |>
as_tibble() |>
mutate(trait = "thallus dry matter content")
#...........Incorporating Species as a Prediction Term...........
tdmc_ir_predictions_species <- ggpredict(tdmc_reg_ir, terms = c("weight_r", "species")) |>
as_tibble() |>
rename(species = group)
# look at the data frame
tdmc_ir_predictions## # A tibble: 29,051 × 7
## x predicted std.error conf.low conf.high group trait
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
## 1 0 -860. 390. -1630. -89.4 1 thallus dry matter content
## 2 1 -858. 390. -1628. -87.6 1 thallus dry matter content
## 3 2 -856. 390. -1626. -85.8 1 thallus dry matter content
## 4 3 -854. 389. -1624. -84.1 1 thallus dry matter content
## 5 4 -852. 389. -1623. -82.3 1 thallus dry matter content
## 6 5 -851. 389. -1621. -80.6 1 thallus dry matter content
## 7 6 -849. 389. -1619. -78.8 1 thallus dry matter content
## 8 7 -847. 389. -1617. -77.1 1 thallus dry matter content
## 9 8 -845. 389. -1615. -75.3 1 thallus dry matter content
## 10 9 -843. 389. -1613. -73.6 1 thallus dry matter content
## # ℹ 29,041 more rowsThickness
#........Predictions for the Thickness Model Without Species........
thick_ir_predictions <- ggpredict(thick_reg_ir_nospecies, terms = "thick_r[0:2]") |>
mutate(trait = "thickness")
# look at the data frame
thick_ir_predictions## # Predicted values of thick_i
##
## : 1
##
## thick_r | Predicted | 95% CI | trait
## --------------------------------------------
## 0 | 0.13 | 0.04, 0.23 | thickness
## 1 | 1.06 | 0.92, 1.19 | thickness
## 2 | 1.98 | 1.68, 2.28 | thickness
##
## Adjusted for:
## * individual_id = 0 (population-level)Volume
#........Predictions for the Volume Model With Species........
volume_ir_predictions <- ggpredict(volume_reg_ir, terms = "vol_r[0:30]") |>
mutate(trait = "volume")
#...........Incorporating Species as a Prediction Term...........
volume_ir_predictions_species <- ggpredict(volume_reg_ir, terms = c("vol_r", "species")) |>
as_tibble() |>
rename(species = group)
# look at the data frame
volume_ir_predictions## # Predicted values of vol_i
##
## vol_r | Predicted | 95% CI | trait
## -----------------------------------------
## 0 | 1.30 | -0.44, 3.05 | volume
## 1 | 2.74 | 1.06, 4.41 | volume
## 2 | 4.17 | 2.56, 5.77 | volume
## 3 | 5.60 | 4.05, 7.15 | volume
## 4 | 7.03 | 5.53, 8.53 | volume
## 5 | 8.46 | 7.00, 9.92 | volume
## 6 | 9.89 | 8.46, 11.33 | volume
##
## Adjusted for:
## * species = BF
## * individual_id = 0 (population-level)##b. Saving predictions for models with herbarium traits as predictor variables
Height
#........Predictions for the Height Model Without Species........
height_ih_predictions <- ggpredict(height_reg_ih_nospecies, terms = "height_h[0:45]")
height_ih_predictions## # Predicted values of height_i
##
## height_h | Predicted | 95% CI
## -----------------------------------
## 0 | 0.26 | -0.66, 1.18
## 6 | 7.08 | 6.43, 7.74
## 11 | 12.77 | 12.28, 13.26
## 17 | 19.59 | 19.14, 20.04
## 23 | 26.42 | 25.80, 27.03
## 28 | 32.10 | 31.28, 32.93
## 34 | 38.93 | 37.82, 40.04
## 45 | 51.44 | 49.77, 53.11
##
## Adjusted for:
## * individual_id = 0 (population-level)Width
#........Predictions for the Width Model Without Species........
width_ih_predictions <- ggpredict(width_reg_ih_nospecies, terms = "width_h[0:30]") #max width = 24
# look at the data frame
width_ih_predictions## # Predicted values of width_i
##
## width_h | Predicted | 95% CI
## ----------------------------------
## 0 | 0.09 | -0.08, 0.26
## 4 | 4.91 | 4.79, 5.04
## 7 | 8.53 | 8.34, 8.72
## 11 | 13.36 | 13.04, 13.68
## 15 | 18.19 | 17.72, 18.65
## 19 | 23.01 | 22.40, 23.62
## 23 | 27.84 | 27.08, 28.59
## 30 | 36.28 | 35.27, 37.30
##
## Adjusted for:
## * individual_id = 0 (population-level) #...........Incorporating Species as a Prediction Term...........
width_ih_predictions_species <- ggpredict(width_reg_ih, terms = c("width_h", "species")) |>
as_tibble() |>
rename(species = group)Perimeter
#........Predictions for the Perimeter Model With Species........
perimeter_ih_predictions <- ggpredict(perimeter_reg_ih, terms = "per_h[0:375]")
#...........Incorporating Species as a Prediction Term...........
perimeter_ih_predictions_species <- ggpredict(perimeter_reg_ih, terms = c("per_h", "species")) |>
as_tibble() |>
rename(species = group)
# look at the data frame
perimeter_ih_predictions_species## # A tibble: 180 × 6
## x predicted std.error conf.low conf.high species
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 0 2.60 13.8 -24.7 29.9 BF
## 2 0 -4.22 14.1 -32.1 23.7 CC
## 3 0 -16.8 19.6 -55.4 21.9 CO
## 4 0 9.24 16.0 -22.3 40.8 CYOS
## 5 0 36.9 11.5 14.3 59.6 DP
## 6 0 0.575 14.0 -27.0 28.2 EGME
## 7 0 3.65 18.2 -32.4 39.8 GR
## 8 0 -17.0 45.8 -108. 73.6 NAND
## 9 0 150. 39.6 71.5 228. POLA
## 10 0 1.15 12.5 -23.6 25.9 R
## # ℹ 170 more rowsSurface Area
#........Predictions for the Surface Area Model With Species........
surf_ih_predictions <- ggpredict(surf_reg_ih, terms = "surf_h[0:430]")
#...........Incorporating Species as a Prediction Term...........
surf_ih_predictions_species <- ggpredict(surf_reg_ih, terms = c("surf_h", "species")) |>
as_tibble() |>
rename(species = group)
# look at the data frame
surf_ih_predictions_species## # A tibble: 150 × 6
## x predicted std.error conf.low conf.high species
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 0 -0.463 7.16 -14.6 13.7 BF
## 2 0 1.23 5.81 -10.3 12.7 CC
## 3 0 1.87 6.10 -10.2 13.9 CO
## 4 0 2.35 7.58 -12.7 17.4 CYOS
## 5 0 14.7 4.49 5.82 23.6 DP
## 6 0 0.889 5.93 -10.8 12.6 EGME
## 7 0 -2.46 12.7 -27.6 22.6 GR
## 8 0 1.21 12.1 -22.7 25.2 NAND
## 9 0 -87.7 24.9 -137. -38.4 POLA
## 10 0 0.239 7.02 -13.6 14.1 R
## # ℹ 140 more rowsWeight
#........Predictions for the Weight Model With Species........
tdmc_ih_predictions <- ggpredict(tdmc_reg_ih, terms = "weight_h[0:29050]")
#...........Incorporating Species as a Prediction Term...........
tdmc_ih_predictions_species <- ggpredict(tdmc_reg_ih, terms = c("weight_h", "species")) |>
as_tibble() |>
rename(species = group)
# look at the data frame
tdmc_ih_predictions## # Predicted values of weight_i
##
## weight_h | Predicted | 95% CI
## -------------------------------------------
## 0 | -112.64 | -1296.85, 1071.57
## 3631 | 28230.75 | 26842.96, 29618.54
## 7263 | 56581.94 | 53600.77, 59563.12
## 10894 | 84925.33 | 80204.85, 89645.81
## 14525 | 113268.72 | 106779.23, 119758.20
## 18156 | 141612.10 | 133342.94, 149881.26
## 21787 | 169955.49 | 159901.65, 180009.33
## 29050 | 226650.07 | 213019.24, 240280.90
##
## Adjusted for:
## * species = BF
## * individual_id = 0 (population-level)tdmc_ih_predictions_species## # A tibble: 170 × 6
## x predicted std.error conf.low conf.high species
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 0 -113. 599. -1297. 1072. BF
## 2 0 93.3 708. -1307. 1493. CC
## 3 0 -136. 743. -1605. 1333. CO
## 4 0 115. 1127. -2114. 2344. CYOS
## 5 0 43.3 528. -1001. 1088. DP
## 6 0 55.5 725. -1378. 1489. EGME
## 7 0 -1659. 1345. -4318. 1000. GR
## 8 0 2199. 1611. -986. 5384. NAND
## 9 0 17.2 1706. -3356. 3391. POLA
## 10 0 111. 691. -1257. 1478. R
## # ℹ 160 more rowsThickness
#........Predictions for the Thickness Model Without Species........
thick_ih_predictions <- ggpredict(thick_reg_ih_nospecies, terms = "thick_h[0:2]")
# look at the data frame
thick_ih_predictions## # Predicted values of thick_i
##
## thick_h | Predicted | 95% CI
## --------------------------------
## 0 | 0.20 | 0.12, 0.29
## 1 | 1.16 | 1.00, 1.32
## 2 | 2.12 | 1.78, 2.47
##
## Adjusted for:
## * individual_id = 0 (population-level) #...........Incorporating Species as a Prediction Term...........
thick_ih_predictions_species <- ggpredict(thick_reg_ih, terms = c("thick_h", "species")) |>
as_tibble() |>
rename(species = group)Volume
#........Predictions for the Volume Model With Species........
volume_ih_predictions <- ggpredict(volume_reg_ih, terms = "vol_h[0:30]")
#...........Incorporating Species as a Prediction Term...........
volume_ih_predictions_species <- ggpredict(volume_reg_ih, terms = c("vol_h", "species")) |>
as_tibble() |>
rename(species = group)
# look at the data frame
volume_ih_predictions## # Predicted values of vol_i
##
## vol_h | Predicted | 95% CI
## ----------------------------------
## 0 | -0.24 | -2.30, 1.82
## 4 | 24.45 | 22.50, 26.39
## 7 | 42.96 | 39.11, 46.81
## 11 | 67.65 | 60.95, 74.34
## 15 | 92.33 | 82.73, 101.93
## 19 | 117.02 | 104.49, 129.54
## 23 | 141.70 | 126.24, 157.16
## 30 | 184.90 | 164.30, 205.51
##
## Adjusted for:
## * species = BF
## * individual_id = 0 (population-level)volume_ih_predictions_species## # A tibble: 209 × 6
## x predicted std.error conf.low conf.high species
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 0.1 0.378 1.01 -1.63 2.38 BF
## 2 0.1 0.362 1.14 -1.89 2.61 CC
## 3 0.1 0.296 0.925 -1.53 2.12 CO
## 4 0.1 1.02 1.87 -2.69 4.72 CYOS
## 5 0.1 4.95 4.40 -3.74 13.6 DL
## 6 0.1 1.10 0.765 -0.410 2.62 DP
## 7 0.1 0.686 0.799 -0.894 2.27 EGME
## 8 0.1 -1.49 2.47 -6.38 3.39 GR
## 9 0.1 1.25 2.02 -2.75 5.25 NAND
## 10 0.1 3.32 3.02 -2.65 9.29 POLA
## # ℹ 199 more rows6. Visualizing Models
This section includes code for scatterplots of model predictions using ggplot. Every model has one scatterplot with initial trait values on the Y axis and treatment (rehydrated/herbarium) trait values on the X axis. Models that incorporate a species interaction have an additional scatterplot faceted by species.
##a. Visualizing models with Rehydrated traits on the X axis
Height
#.....Scatterplot of height_r(x) as Predictor of height_i(y).....
height_ir_plot <- ggplot() +
# plot the raw data as points
geom_point(data = height_all,
aes(x = height_r,
y = height_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = height_ir_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = height_ir_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Rehydrate height (mm)",
y = "Initial height (mm)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 50, by = 10),
limits = c(0, 50)) +
scale_y_continuous(breaks = seq(from = 0, to = 50, by = 10),
limits = c(0, 50))+
#add equation + r squared
annotate("text", x = 10, y = 45, label= "y = 1.02x + 0.0407
R squared = 0.996")
height_ir_plot## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
### Width
#.....Scatterplot of width_r(x) as Predictor of width_i(y).....
width_ir_plot <- ggplot() +
# plot the raw data as points
geom_point(data = width_all,
aes(x = width_r,
y = width_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = width_ir_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = width_ir_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Rehydrate Width (mm)",
y = "Initial width (mm)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 30, by = 5),
limits = c(0, 30)) +
scale_y_continuous(breaks = seq(from = 0, to = 30, by = 5),
limits = c(0, 30))+
#add equation + r squared
annotate("text", x = 10, y = 25, label= "y = 0.994x + 0.0865
R squared = 0.996")
width_ir_plot
### Thickness
#.....Scatterplot of thickness_r(x) as Predictor of thickness_i(y).....
thick_ir_plot <- ggplot() +
# plot the raw data as points
geom_point(data = thickness_all,
aes(x = thick_r,
y = thick_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = thick_ir_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = thick_ir_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Rehydrate thickness (mm)",
y = "Initial thickness (mm)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 2, by = 0.25),
limits = c(0, 2)) +
scale_y_continuous(breaks = seq(from = 0, to = 2, by = 0.25),
limits = c(0, 2))+
#add equation + r squared
annotate("text", x = 0.75, y = 1.75, label= "y = 0.922x + 0.134
R squared = 0.882")
thick_ir_plot## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).
Perimeter
#.....Scatterplot of perimeter_r(x) as Predictor of perimeter_i(y).....
perimeter_ir_plot <- ggplot() +
# plot the raw data as points
geom_point(data = perimeter_all,
aes(x = per_r,
y = per_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = perimeter_ir_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = perimeter_ir_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Rehydrate perimeter (mm)",
y = "Initial perimeter (mm)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 375, by = 50),
limits = c(0, 375)) +
scale_y_continuous(breaks = seq(from = 0, to = 375, by = 50),
limits = c(0, 375))+
#add equation + r squared
annotate("text", x = 75, y = 300, label= "y = 1.02x + 0.0407
R squared = 0.996")
perimeter_ir_plot## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
#.Scatterplot of height_r(x) as Predictor of height_i(y) with a species facet.
perimeter_ir_plot_species <- ggplot() +
geom_point(data = perimeter_all,
aes(x = per_r,
y = per_i,
color = species),
alpha = 0.5,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence intervals
geom_ribbon(data = perimeter_ir_predictions_species,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high,
fill = species,
group = species),
alpha = 0.2) +
# plot the model predictions
geom_line(data = perimeter_ir_predictions_species,
aes(x = x,
y = predicted,
color = species,
group = species)) +
# labels
labs(x = "Rehydrate perimeter (mm)",
y = "Initial perimeter (mm)") +
# theme
theme(legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(size = 12)) +
# facet
facet_wrap(~ species, nrow = 2,
scales = "free")
perimeter_ir_plot_species## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
### Surface Area
#.....Scatterplot of surf_r(x) as Predictor of surf_i(y).....
surf_ir_plot <- ggplot() +
# plot the raw data as points
geom_point(data = surf_all,
aes(x = surf_r,
y = surf_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = surf_ir_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = surf_ir_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Rehydrate Surface Area (mm^2)",
y = "Initial Surface Area (mm^2)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 450, by = 75),
limits = c(0, 450)) +
scale_y_continuous(breaks = seq(from = 0, to = 450, by = 75),
limits = c(0, 450))
surf_ir_plot## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
#.Scatterplot of surf_r(x) as Predictor of surf_i(y) with a species facet.
surf_ir_plot_species <- ggplot() +
geom_point(data = surf_all,
aes(x = surf_r,
y = surf_i,
color = species),
alpha = 0.5,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence intervals
geom_ribbon(data = surf_ir_predictions_species,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high,
fill = species,
group = species),
alpha = 0.2) +
# plot the model predictions
geom_line(data = surf_ir_predictions_species,
aes(x = x,
y = predicted,
color = species,
group = species)) +
# labels
labs(x = "Rehydrate surface area (mm^2)",
y = "Initial surface area (mm^2)") +
# theme
theme(legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(size = 12)) +
# facet
facet_wrap(~ species, nrow = 2,
scales = "free")
surf_ir_plot_species
### Weight
#.....Scatterplot of weight_r(x) as Predictor of weight_i(y).....
tdmc_ir_plot <- ggplot() +
# plot the raw data as points
geom_point(data = tdmc_all,
aes(x = weight_r,
y = weight_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = tdmc_ir_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = tdmc_ir_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Rehydrate weight (g)",
y = "Initial weight (g)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 4850, by = 810),
limits = c(0, 4850)) +
scale_y_continuous(breaks = seq(from = 0, to = 4850, by = 810),
limits = c(0, 4850))
tdmc_ir_plot## Warning: Removed 54 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning: Removed 26370 rows containing missing values or values outside the scale range
## (`geom_line()`).
#.Scatterplot of weight_r(x) as Predictor of weight_i(y) with a species facet.
tdmc_ir_plot_species <- ggplot() +
geom_point(data = tdmc_all,
aes(x = weight_r,
y = weight_i,
color = species),
alpha = 0.5,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence intervals
geom_ribbon(data = tdmc_ir_predictions_species,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high,
fill = species,
group = species),
alpha = 0.2) +
# plot the model predictions
geom_line(data = tdmc_ir_predictions_species,
aes(x = x,
y = predicted,
color = species,
group = species)) +
# labels
labs(x = "Rehydrate weight (g)",
y = "Initial weight (g)") +
# theme
theme(legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(size = 12)) +
# facet
facet_wrap(~ species, nrow = 2,
scales = "free")
tdmc_ir_plot_species## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
### Volume
#.....Scatterplot of vol_r(x) as Predictor of vol_i(y).....
volume_ir_plot <- ggplot() +
# plot the raw data as points
geom_point(data = volume_all,
aes(x = vol_r,
y = vol_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = volume_ir_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = volume_ir_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Rehydrate volume (mL)",
y = "Initial volume (mL)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 15, by = 2.5),
limits = c(0, 15)) +
scale_y_continuous(breaks = seq(from = 0, to = 15, by = 2.5),
limits = c(0, 15))
volume_ir_plot## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
#.Scatterplot of volume_r(x) as Predictor of volume_i(y) with a species facet.
volume_ir_plot_species <- ggplot() +
geom_point(data = volume_all,
aes(x = vol_r,
y = vol_i,
color = species),
alpha = 0.5,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence intervals
geom_ribbon(data = volume_ir_predictions_species,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high,
fill = species,
group = species),
alpha = 0.2) +
# plot the model predictions
geom_line(data = volume_ir_predictions_species,
aes(x = x,
y = predicted,
color = species,
group = species)) +
# labels
labs(x = "Rehydrate volume (mL)",
y = "Initial volume (mL)") +
# theme
theme(legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(size = 12)) +
# facet
facet_wrap(~ species, nrow = 2,
scales = "free")
volume_ir_plot_species
##b. Visualizing models with herbarium traits on the X axis
Height
#....Scatterplot of height_h(x) as a Predictor of height_i(y)....
height_ih_plot <- ggplot() +
# plot the raw data as points
geom_point(data = height_all,
aes(x = height_h,
y = height_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = height_ih_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = height_ih_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Herbarium height (mm)",
y = "Initial height (mm)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 50, by = 10),
limits = c(0, 50)) +
scale_y_continuous(breaks = seq(from = 0, to = 50, by = 10),
limits = c(0, 50))+
#add equation + r squared
annotate("text", x = 10, y = 45, label= "y = 1.14x + 0.258
R squared = 0.929")
height_ih_plot## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
Width
#....Scatterplot of width_h(x) as a Predictor of width_i(y)....
width_ih_plot <- ggplot() +
# plot the raw data as points
geom_point(data = width_all,
aes(x = width_h,
y = width_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = width_ih_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = width_ih_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Herbarium width (mm)",
y = "Initial width (mm)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 20, by = 5),
limits = c(0, 20)) +
scale_y_continuous(breaks = seq(from = 0, to = 20, by = 5),
limits = c(0, 20))
width_ih_plot## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_line()`).
#.Scatterplot of width_h(x) as Predictor of width_i(y) with a species facet.
width_ih_plot_species <- ggplot() +
geom_point(data = width_all,
aes(x = width_h,
y = width_i,
color = species),
alpha = 0.5,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence intervals
geom_ribbon(data = width_ih_predictions_species,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high,
fill = species,
group = species),
alpha = 0.2) +
# plot the model predictions
geom_line(data = width_ih_predictions_species,
aes(x = x,
y = predicted,
color = species,
group = species)) +
# labels
labs(x = "Herbarium width (mm)",
y = "Initial width (mm)") +
# theme
theme(legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(size = 12)) +
# facet
facet_wrap(~ species, nrow = 2,
scales = "free")
width_ih_plot_species
### Thickness
#....Scatterplot of thick_h(x) as a Predictor of thick_i(y)....
thick_ih_plot <- ggplot() +
# plot the raw data as points
geom_point(data = thickness_all,
aes(x = thick_h,
y = thick_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = thick_ih_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = thick_ih_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Rehydrate thickness (mm)",
y = "Initial thickness (mm)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 2, by = 0.25),
limits = c(0, 2)) +
scale_y_continuous(breaks = seq(from = 0, to = 2, by = 0.25),
limits = c(0, 2))
thick_ih_plot## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
#.Scatterplot of thick_h(x) as Predictor of thick_i(y) with a species facet.
thick_ih_plot_species <- ggplot() +
geom_point(data = thickness_all,
aes(x = thick_h,
y = thick_i,
color = species),
alpha = 0.5,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence intervals
geom_ribbon(data = thick_ih_predictions_species,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high,
fill = species,
group = species),
alpha = 0.2) +
# plot the model predictions
geom_line(data = thick_ih_predictions_species,
aes(x = x,
y = predicted,
color = species,
group = species)) +
# labels
labs(x = "Herbarium thickness (mm)",
y = "Initial thickness (mm)") +
# theme
theme(legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(size = 12)) +
# facet
facet_wrap(~ species, nrow = 2,
scales = "free")
thick_ih_plot_species
### Perimeter
#....Scatterplot of per_h(x) as a Predictor of per_i(y)....
perimeter_ih_plot <- ggplot() +
# plot the raw data as points
geom_point(data = perimeter_all,
aes(x = per_h,
y = per_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = perimeter_ih_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = perimeter_ih_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Herbarium perimeter (mm)",
y = "Initial perimeter (mm)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 375, by = 50),
limits = c(0, 375)) +
scale_y_continuous(breaks = seq(from = 0, to = 375, by = 50),
limits = c(0, 375))
perimeter_ih_plot## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning: Removed 63 rows containing missing values or values outside the scale range
## (`geom_line()`).
#.Scatterplot of per_h(x) as Predictor of per_i(y) with a species facet.
perimeter_ih_plot_species <- ggplot() +
geom_point(data = perimeter_all,
aes(x = per_h,
y = per_i,
color = species),
alpha = 0.5,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence intervals
geom_ribbon(data = perimeter_ih_predictions_species,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high,
fill = species,
group = species),
alpha = 0.2) +
# plot the model predictions
geom_line(data = perimeter_ih_predictions_species,
aes(x = x,
y = predicted,
color = species,
group = species)) +
# labels
labs(x = "Herbarium perimeter (mm)",
y = "Initial perimeter (mm)") +
# theme
theme(legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(size = 12)) +
# facet
facet_wrap(~ species, nrow = 2,
scales = "free")
perimeter_ih_plot_species## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
### Surface Area
#....Scatterplot of surf_h(x) as a Predictor of surf_i(y)....
surf_ih_plot <- ggplot() +
# plot the raw data as points
geom_point(data = surf_all,
aes(x = surf_h,
y = surf_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = surf_ih_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = surf_ih_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Herbarium Surface Area (mm^2)",
y = "Initial Surface Area (mm^2)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 450, by = 75),
limits = c(0, 450)) +
scale_y_continuous(breaks = seq(from = 0, to = 450, by = 75),
limits = c(0, 450))
surf_ih_plot## Warning: Removed 159 rows containing missing values or values outside the scale range
## (`geom_line()`).
#.Scatterplot of surf_h(x) as Predictor of surf_i(y) with a species facet.
surf_ih_plot_species <- ggplot() +
geom_point(data = surf_all,
aes(x = surf_r,
y = surf_i,
color = species),
alpha = 0.5,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence intervals
geom_ribbon(data = surf_ih_predictions_species,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high,
fill = species,
group = species),
alpha = 0.2) +
# plot the model predictions
geom_line(data = surf_ih_predictions_species,
aes(x = x,
y = predicted,
color = species,
group = species)) +
# labels
labs(x = "Herbarium Surface Area (mm^2)",
y = "Initial Surface Area (mm^2)") +
# theme
theme(legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(size = 12)) +
# facet
facet_wrap(~ species, nrow = 2,
scales = "free")
surf_ih_plot_species
### Weight
#....Scatterplot of weight_h(x) as a Predictor of weight_i(y)....
tdmc_ih_plot <- ggplot() +
# plot the raw data as points
geom_point(data = tdmc_all,
aes(x = weight_h,
y = weight_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = tdmc_ih_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = tdmc_ih_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Herbarium weight (g)",
y = "Initial weight (g)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 4850, by = 810),
limits = c(0, 4850)) +
scale_y_continuous(breaks = seq(from = 0, to = 4850, by = 810),
limits = c(0, 4850))
tdmc_ih_plot## Warning: Removed 54 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning: Removed 28430 rows containing missing values or values outside the scale range
## (`geom_line()`).
#.Scatterplot of weight_h(x) as Predictor of weight_i(y) with a species facet.
tdmc_ih_plot_species <- ggplot() +
geom_point(data = tdmc_all,
aes(x = weight_h,
y = weight_i,
color = species),
alpha = 0.5,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence intervals
geom_ribbon(data = tdmc_ih_predictions_species,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high,
fill = species,
group = species),
alpha = 0.2) +
# plot the model predictions
geom_line(data = tdmc_ih_predictions_species,
aes(x = x,
y = predicted,
color = species,
group = species)) +
# labels
labs(x = "Herbarium weight (g)",
y = "Initial weight (g)") +
# theme
theme(legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(size = 12)) +
# facet
facet_wrap(~ species, nrow = 2,
scales = "free")
tdmc_ih_plot_species## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
### Volume
#....Scatterplot of vol_h(x) as a Predictor of vol_i(y)....
volume_ih_plot <- ggplot() +
# plot the raw data as points
geom_point(data = volume_all,
aes(x = vol_h,
y = vol_i),
alpha = 0.5,
color = model_col,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence interval
geom_ribbon(data = volume_ih_predictions,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high),
alpha = 0.2,
fill = model_col) +
# plot the prediction line
geom_line(data = volume_ih_predictions,
aes(x = x,
y = predicted),
color = model_col,
linewidth = 1) +
# labels
labs(x = "Herbarium volume (mL)",
y = "Initial volume (mL)") +
# controlling axes to make the plot look square
# this makes it easier to see the difference between the model and the 1:1 line
scale_x_continuous(breaks = seq(from = 0, to = 15, by = 2.5),
limits = c(0, 15)) +
scale_y_continuous(breaks = seq(from = 0, to = 15, by = 2.5),
limits = c(0, 15))
volume_ih_plot## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning: Removed 29 rows containing missing values or values outside the scale range
## (`geom_line()`).
#.Scatterplot of vol_h(x) as Predictor of vol_i(y) with a species facet.
volume_ih_plot_species <- ggplot() +
geom_point(data = volume_all,
aes(x = vol_h,
y = vol_i,
color = species),
alpha = 0.5,
shape = 21) +
# plot a 1:1 reference line
geom_abline(slope = 1,
intercept = 0,
linetype = 2,
linewidth = 1,
color = ref_line_col) +
# plot the confidence intervals
geom_ribbon(data = volume_ih_predictions_species,
aes(x = x,
y = predicted,
ymin = conf.low,
ymax = conf.high,
fill = species,
group = species),
alpha = 0.2) +
# plot the model predictions
geom_line(data = volume_ih_predictions_species,
aes(x = x,
y = predicted,
color = species,
group = species)) +
# labels
labs(x = "Herbarium volume (mL)",
y = "Initial volume (mL)") +
# theme
theme(legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(size = 12)) +
# facet
facet_wrap(~ species, nrow = 2,
scales = "free")
volume_ih_plot_species
c. Visualizing summary graphs for rehydrate models
This section includes visualizations that plot multiple trait models on a single graph
Scatterplot faceted by trait
This scatterplot visualizes predictions for height, width, perimeter, weight, surface area, thickness, and volume as a function of inital trait values, faceted by trait.
# surf_ir_predictions, perimeter_ir_predictions, tdmc_ir_predictions, height_ir_predictions, width_ir_predictions, thick_ir_predictions, volume_ir_predictions
#.........................Data Wrangling.........................
#creating a list of all the predictions for rehydrate traits
init_rehy_df <- rbind(surf_ir_predictions, # rbind stacks data frames on top of each other if they have the same columns
perimeter_ir_predictions,
tdmc_ir_predictions,
height_ir_predictions,
width_ir_predictions,
thick_ir_predictions,
volume_ir_predictions) |>
drop_na() # not sure why there are NA rows at the end of the data frame, but they come from one of the prediction data frames
#.........A scatterplot with all seven model predictions.........
ggplot(data = init_rehy_df, # using the data frame created above
aes(x = x, # x-axis
y = predicted )) + # y-axis
geom_abline(slope = 1, # adding a dashed grey 1:1 reference line
intercept = 0,
color = "grey",
linetype = 2) +
# geom_ribbon(aes(fill = trait, # adding 95% CI
# ymin = conf.low,
# ymax = conf.high),
# alpha = 0.2) +
geom_line(aes(color = trait),
linewidth = 1) + # plotting model predictions
labs(x = "Rehydrated value", # renaming the axes
y = "Initial value") +
theme(legend.position = "none") +
facet_wrap(~trait, scales = "free")
Highlighting 1:1 relationships
This scatterplot visualizes predictions for height, width, perimeter, surface area, thickness, and volume as a function of inital trait values. Prediction lines are on the same graph, and models that follow a 1:1 relationship are highlighted.
#..............................Setup.............................
#Using the dataframe with all our rehydrate trait predictions to assign an opacity
#value to all traits. 3 signifies a close 1:1 relationship while 0.1 is a weak 1:1.
init <- init_rehy_df
init23 <- init[init$trait != "thallus dry matter content"]
init23$opacity <- ifelse(init23$trait == "height", 3,
ifelse(init23$trait == "width", 3,
ifelse(init23$trait == "perimeter", 0.1,
ifelse(init23$trait == "surface area", 0.1,
ifelse(init23$trait == "thickness", 0.1,
ifelse(init23$trait == "volume", 0.1, NA ))))))
#........manually offsetting height and width line labels........
label_df <- init23 |>
filter(x == max(x)) |>
mutate(label = stringr::str_to_title(trait),
nudge_y = case_when(
trait == "height" ~ 0.7,
trait == "width" ~ -0.3,
TRUE ~ 0
))
#........................generating titles.......................
title <- "Two predictive trait models based on rehydrated values closely follow a 1:1 relationship"
subtitle <- "The traits that follow this pattern are Height and Width"
#..............assigning each trait a unique color...............
showtext_auto()
trait_palette <- c("surface area" = "#ffbf00",
"perimeter" = "#be8333",
# "thallus dry matter content" = "#54662c",
"height" = "#009bb0",
"width" = "#124c54",
"thickness" = "#368000",
"volume" = "#b4450e")
#........scatterplot with all rehydrate trait predictions........
#......rehydrated trait value (x) & initial trait values (y).....
#..excludes tdmc, and highlights traits with 1:1 relationships...
plot_ir_1s <- init23 |>
ggplot(aes(x, predicted, color = trait, alpha = opacity))+
geom_line()+
geom_abline(slope = 1, # adding a dashed grey 1:1 reference line
intercept = 0,
color = "grey",
linetype = 2) +
# geom_ribbon(aes(fill = trait, # adding 95% CI (uncomment this if you want to see it)
# ymin = conf.low,
# ymax = conf.high),
# alpha = 0.2) +
geom_line(
linewidth = 0.5) + # plotting model predictions
labs(x = "Rehydrated Value", # renaming the axes
y = "Predicted Initial Value",
title = title,
subtitle = subtitle) +
scale_color_manual(values = trait_palette)+ #assigning values a color based on trait from our palette
theme(legend.position = "none",
plot.background = element_blank())+
#formatting line labels
geom_text(
data = init23 |> filter(x == max(x)) |> mutate(label = stringr::str_to_title(trait)),
aes(label = label, alpha = opacity),
hjust = 0,
nudge_x = 0.1,
nudge_y = label_df$nudge_y,
size = 5,
family = "lex",
#alpha = opacity
#fontface = "light"
)+
#formatting text and margins
coord_cartesian(xlim = c(0, 7.2),
clip = "off")+
theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 1, "cm"))+
theme(axis.title.x = element_text(family = "lex",
size = 14,
hjust = 0.5,
margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")),
axis.title.y = element_text(family = "lex",
size = 14,
hjust = 0.5,
margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")),
plot.title = element_text(family = "lex",
size = 17,
hjust = 0.5,
margin = margin(t = 0, r = 1, b = 0.3, l = 0, "cm")),
plot.subtitle = element_text(family = "lex",
size = 15,
hjust = 0.5,
margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")))+
scale_x_continuous(expand = c(0, NA))+
scale_y_continuous(limits = c(0.8,NA))
plot_ir_1s## Warning in y + params$y: longer object length is not a multiple of shorter
## object length## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_text()`).
Highlighting relationships that deviate from 1:1
This includes code for a scatterplot of height, width, perimeter, surface area, thickness, and volume on the same graph. Models that deviate from a 1:1 relationship are highlighted.
#.........................assigning opacity values (manually).........................
#traits that follow a 1:1 relationship are assigned 1, while deviants from 1:1
#are assigned 3.
showtext_auto()
init <- init_rehy_df
init2 <- as.data.frame(init) |>
mutate(trait = as.factor(trait))
init2 <- init[init$trait != "thallus dry matter content"]
init2$opacity <- ifelse(init2$trait == "height", 1,
ifelse(init2$trait == "width", 1,
ifelse(init2$trait == "perimeter", 3,
ifelse(init2$trait == "surface area", 3,
ifelse(init2$trait == "thickness", 1,
ifelse(init2$trait == "volume", 3, NA ))))))
#.........................titles and palettes.........................
title <- "Predictive models based on rehydrated measurements vary by trait"
subtitle <- "The traits that deviate the most are Perimeter, Volume, and Surface Area"
showtext_auto()
trait_palette <- c("surface area" = "#ffbf00",
"perimeter" = "#be8333",
# "thallus dry matter content" = "#54662c",
"height" = "#009bb0",
"width" = "#124c54",
"thickness" = "#368000",
"volume" = "#b4450e")
#.........................plotting.........................
#......rehydrate trait values (x) & initial trait values (y).....
plot_ir_deviants <- init2 |>
ggplot(aes(x, predicted, color = trait, alpha = opacity))+
geom_line()+
geom_abline(slope = 1, # adding a dashed grey 1:1 reference line
intercept = 0,
color = "grey",
linetype = 2) +
# geom_ribbon(aes(fill = trait, # adding 95% CI (uncomment this if you want to see it)
# ymin = conf.low,
# ymax = conf.high),
# alpha = 0.2) +
geom_line(
linewidth = 1.5) + # plotting model predictions
labs(x = "Rehydrated Value", # renaming the axes
y = "Predicted Initial Value",
title = title,
subtitle = subtitle) +
scale_color_manual(values = trait_palette)+ #assigning palette colors to traits
theme(legend.position = "none",
plot.background = element_blank())+
geom_text( #formatting model labels
data = init2 |> filter(x == max(x)) |> mutate(label = stringr::str_to_title(trait)),
aes(label = label),
hjust = 0,
nudge_x = 0.1,
size = 5,
family = "lex"
#fontface = "light"
)+
#manually adjusting text elements and plot dimensions
coord_cartesian(xlim = c(0, 7.2),
clip = "off")+
theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 1, "cm"))+
theme(axis.title.x = element_text(family = "lex",
size = 14,
hjust = 0.5,
margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")),
axis.title.y = element_text(family = "lex",
size = 14,
hjust = 0.5,
margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")),
plot.title = element_text(family = "lex",
size = 17,
hjust = 0.5,
margin = margin(t = 0, r = 0, b = 0.3, l = 0, "cm")),
plot.subtitle = element_text(family = "lex",
size = 15,
hjust = 0.5,
margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")))+
scale_x_continuous(expand = c(0, NA))+
scale_y_continuous(limits = c(0.8,NA))
plot_ir_deviants## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_text()`).